rm(list = ls())

## SET YOU FILE PATH
setwd("~/Dropbox/Immigration/Analysis/Policy Diffusion/Replication Files")

library(ggplot2)
library(ggpubr)
library(dplyr)
library(survminer)
library(tidyr)
library(standardize)
library(survival)
library(coxme)
library(frailtySurv)
library(frailtyEM)
library(frailtypack)
library(Hmisc)
library(haven)
library(foreign)
library(coefplot)
library(stargazer)
library(simPH)
library(coefplot)
library(xtable)
library(sp)
library(sf)
require(rgdal)
library(tmap)
library(tmaptools)
library(sp)
library(leaflet)
library(haven)
library(ggplot2)
library(shinyjs)
library(RColorBrewer)
library(openssl)
library(texreg)

data <- read_dta("dwrap_index.dta")

data.africa <- subset(data, africa == 1)
data.1980 <- subset(data, year >= 1980)
data.indp5 <- subset(data, indp5 == 1)
data.law <- subset(data, numdomlaws >= 1)

set.seed(8675309)

## Set Convergence Criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)

coxme.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000)

## TABLE 2

res1 <- matrix(NA,8,8)

m1c <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
               majorannintra_contig_lag1 + 
               
               strata(numchg1)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1c)
res1[1:2,1] <- c(m1c$coef[1],sqrt(m1c$var[1,1]))
res1[3,1] <- summary(m1c)$coefficients[1, 6]
res1[5,1] <- m1c$nevent
res1[6,1] <- m1c$n
res1[7,1] <- m1c$history[[1]]$c.loglik
res1[8,1] <- AIC(m1c)

m2c <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
               majorannintra_contig_lag1 +
               exc1inc2_reg_lag1+
               ihs_aidgdp_5ma_wins+
               polyarchy_lag1 +
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               strata(numchg1)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2c)
res1[1:2,2] <- c(m2c$coef[1],sqrt(m2c$var[1,1]))
res1[3,2] <- summary(m2c)$coefficients[1, 6]
res1[5,2] <- m2c$nevent
res1[6,2] <- m2c$n
res1[7,2] <- m2c$history[[1]]$c.loglik
res1[8,2] <- AIC(m2c)

m3c <- coxph(formula = Surv(ewchg1_start, ewchg1_end, ew_change1) ~
               majorannintra_contig_lag1 + 
               
               strata(ewnumchg1)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3c)
res1[1:2,3] <- c(m3c$coef[1],sqrt(m3c$var[1,1]))
res1[3,3] <- summary(m3c)$coefficients[1, 6]
res1[5,3] <- m3c$nevent
res1[6,3] <- m3c$n
res1[7,3] <- m3c$history[[1]]$c.loglik
res1[8,3] <- AIC(m3c)

m4c <- coxph(formula = Surv(chg5_start, chg5_end, change_50) ~
               majorannintra_contig_lag1 + 
               
               strata(numchg5)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m4c)
res1[1:2,4] <- c(m4c$coef[1],sqrt(m4c$var[1,1]))
res1[3,4] <- summary(m4c)$coefficients[1, 6]
res1[5,4] <- m4c$nevent
res1[6,4] <- m4c$n
res1[7,4] <- m4c$history[[1]]$c.loglik
res1[8,4] <- AIC(m4c)

m5c <- coxph(formula = Surv(ewchg5_start, ewchg5_end, ew_change_5) ~
               majorannintra_contig_lag1 + 
               
               strata(ewnumchg5)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m5c)
res1[1:2,5] <- c(m5c$coef[1],sqrt(m5c$var[1,1]))
res1[3,5] <- summary(m5c)$coefficients[1, 6]
res1[5,5] <- m5c$nevent
res1[6,5] <- m5c$n
res1[7,5] <- m5c$history[[1]]$c.loglik
res1[8,5] <- AIC(m5c)

m6c <- coxph(formula = Surv(chg15_start, chg15_end, change_15) ~
               majorannintra_contig_lag1 + 
               
               strata(numchg15)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m6c)
res1[1:2,6] <- c(m6c$coef[1],sqrt(m6c$var[1,1]))
res1[3,6] <- summary(m6c)$coefficients[1, 6]
res1[5,6] <- m6c$nevent
res1[6,6] <- m6c$n
res1[7,6] <- m6c$history[[1]]$c.loglik
res1[8,6] <- AIC(m6c)

m7c <- coxph(formula = Surv(ewchg15_start, ewchg15_end, ew_change_15) ~
               majorannintra_contig_lag1 + 
               
               strata(ewnumchg15)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m7c)
res1[1:2,7] <- c(m7c$coef[1],sqrt(m7c$var[1,1]))
res1[3,7] <- summary(m7c)$coefficients[1, 6]
res1[5,7] <- m7c$nevent
res1[6,7] <- m7c$n
res1[7,7] <- m7c$history[[1]]$c.loglik
res1[8,7] <- AIC(m7c)

m8c <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
               majorannintra_reg1500_lag1 + 
               
               strata(numchg1)+
               cluster(ccode2)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m8c)
res1[1:2,8] <- c(m8c$coef[1],sqrt(m8c$var[1,1]))
res1[3,8] <- summary(m8c)$coefficients[1, 6]
res1[5,8] <- m8c$nevent
res1[6,8] <- m8c$n
res1[7,8] <- m8c$history[[1]]$c.loglik
res1[8,8] <- AIC(m8c)

res1 <- round(res1,3)

latex(res1, file = "")

## FIGURE A.26

data$conflict <- NULL
data$conflict <- data$majorannintra_contig_lag1

multi1 <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
                  conflict + 
                 
                 strata(numchg1)+
                 cluster(ccode2)+
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(multi1)

data$conflict <- NULL
data$conflict <- data$intrastate_contig_lag1

multi2 <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
                  conflict + 
                  
                  strata(numchg1)+
                  cluster(ccode2)+
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(multi2)

data$conflict <- NULL
data$conflict <- data$minorannintra_contig_lag1

multi3 <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
                  conflict + 
                  
                  strata(numchg1)+
                  cluster(ccode2)+
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(multi3)

data$conflict <- NULL
data$conflict <- data$minorcumintra_contig_lag1

multi4 <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
                  conflict + 
                  
                  strata(numchg1)+
                  cluster(ccode2)+
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(multi4)

data$conflict <- NULL
data$conflict <- data$majorcumintra_contig_lag1

multi5 <- coxph(formula = Surv(chg1_start, chg1_end, change1) ~
                  conflict + 
                  
                  strata(numchg1)+
                  cluster(ccode2)+
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(multi5)

cwplot <- multiplot(multi1, multi2, multi3, multi4, multi5, secret.weapon=T, innerCI = 2, coefficients = c("conflict"), sort = c("magnitude"), names=c("Intense, Recent Conflict \n (1000+ Battle Deaths in Prior Year)", "All Civil Wars", "Low-Intensity, Recent Conflict \n (25-999 Battle Deaths in Prior Year)" , "Low-Intensity, Cumulative Conflict \n (25-999 Total Battle Deaths)", "Intense, Cumulative Conflict \n (1000+ Total Battle Deaths)"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient", ylab="")
cwplot

data$conflict <- NULL
